# bootstrap.auk
# Tim Menzies
# Sept 12, 2012

@include "lib.awk"

function testStatistic(y,z,   m,n,ymu,zmu,tmp1,tmp2,s1,s2,i) {
    m= length(y)
    n= length(z)
    ymu= mean(y)
    zmu= mean(z)
    for(i in y) tmp1 += (y[i] - ymu)^2
    for(i in z) tmp2 += (z[i] - zmu)^2
    s1= tmp1/(m-1)
    s2= tmp2/(n-1)
    return abs((zmu - ymu)/sqrt(s1/m + s2/n))
}
function bootstrap(y,z,b,conf, \
                   tobs,x,xmu,ymu,zmu,i,yhat,zhat,bigger,asl) {
    b = b ? b : 1000
    conf = conf ? conf : 0.05
    tobs = testStatistic(y,z)
    concat(y,z,x)
    xmu = mean(x)
    ymu = mean(y)
    zmu = mean(z)
    for(i in y) yhat[i] = y[i] - ymu + xmu
    for(i in z) zhat[i] = z[i] - zmu + xmu
    for(i=1;i<=b;i++) 
      bigger += bootstrap1(y,yhat,z,zhat) > tobs 
    asl= bigger/b
    return asl < conf
}
function bootstrap1(y,yhat,z,zhat,     ystar,zstar,i) {
    for(i in y) ystar[i] = any(yhat)
    for(i in z) zstar[i] = any(zhat)
    return testStatistic(ystar,zstar) 
}
function _bootstrap(  b,m,m1,s1,n,m2,s2, one,two,i,c,o) {
  args("-b,1000,-m,100,-n,100,-s1,1,-s2,1,-m1,10,-m2,11,-c,0.05",o)
  b  = o["-b"] ;  m = o["-m"];  n  = o["-n"]; c=o["-c"]
  s1 = o["-s1"]; s2 = o["-s2"]; 
  m1 = o["-m1"]; m2 = o["-m2"]
  for(i = 1;i<=m;i++) 
    one[i] = normal(m1,s1)
  for(i = 1;i<=n;i++) 
    two[i] = normal(m2,s2)
  print "["bootstrap(one,two,b,c)"]"
}
